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Abstract 


GRACE data indicate large seasonal variations in gravity that have been shown to be 
to be related to climate-driven fluxes of surface water. Seasonal redistribution of sur- 
face mass deforms the Earth, and our previous study using GRACE data demonstrate 
that annual radial deformations of ±13 mm in the region of Amazon River Basin were 
observed by both GRACE and ten GPS sites in the region. For the GRACE determi- 
nations, we estimate in a least-squares solution for each Stokes coefficient parameters 
that represent the amplitudes of the annual variation. We then filter these parameters 
based on a statistical test that uses the scatter of the postfit residuals. We demonstrate 
by comparison to the GPS amplitudes that this method is more accurate, for this region, 
than Gaussian smoothing. Our model for the temporal behavior of the gravity coeffi- 
cients includes a rate term, and although the time series are noisy, the glacial isostatic 
adjustment signal over Hudson’s Bay can be observed. 
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1 Introduction 


The GRACE mission is providing a unique opportunity for studying the Earth's global gravity vari- 
ations. At seasonal timescales, the most significant mass motions are associated with climate-driven 
transport of water on the surface of the Earth. This redistribution of water on the surface of the Earth 
will act as a change in the surface load and will deform the Earth. 

GRACE is therefore a potential tool for inferring loading-induced deformation of the solid Earth 
over timescales greater than about one month, the approximate time required to produce a global 
gravity field. This temporal coverage may be compared to surveying with the Global Positioning 
System (GPS) in “continuous" mode, in which estimates of site position may be made by combining 
data acquired every 24 hours or less. GRACE recov ers the gravity variations, and the surface mass 
redistributions derived from them, with global coverage of high spatial resolution (typically up to 
degree l = 120). GPS observations are acquired from a generally sparse (but in some places very 
dense) network that is highly heterogeneous, especially in regard to the southern hemisphere and the 
oceans. 

GRACE and GPS are thus complementary approaches for study ing loading-induced deformation 
of the solid Earth. GRACE has no information at £ — 1. for example, whereas GPS has been used for 
several years to detect the deformation at this degree (Blewitt et al., 2001). Recovering deformation 
at higher degrees poses difficulty for GPS. however (Blewitt et al., 2003), whereas for GRACE it is 
natural. 

In Davis et al. (2004) we compared estimates from GRACE and GPS of annual amplitudes of 
radial deformation for South America, and demonstrated that there is a significant correlation. In 
that work, we briefly discussed a statistical approach for filtering the estimated annual amplitudes. 
Below we extend our calculations for the existing GRACE data set, and investigate the accuracy of 
the statistical filtering approach by comparing the results to those obtained using a more standard 
Gaussian smoothing approach. 

2 Calculation of radial deformation 

A GRACE Level-2 data set (Bettadpur, 2003) consists of Stokes coefficients representing an expan- 
sion in spherical harmonics of the Earth’s gravity potential. Each data set results from the reduction 
of approximately one month of GRACE observations consisting primarily of the range between 
the co-orbiting satellites. Corrections are applied for oceanic and atmospheric mass motions on 
timescales shorter than one month (Flechtner, 2003). Changes in the Stokes coefficients can be used 
to calculate changes in the geoid AN((f>, A) at a point on the surface of the Earth with latitude and 
longitude A using 
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where a is the radius of the Earth, P <m (sin <p) are the fully normalized associated Legendre functions 
(Bettadpur, 2003), A Ce m and A S( m are changes in the dimensionless Stokes coefficients at degree 
£ and order m. (The summation begins at £ = 2 because, as mentioned above, GRACE has no 
information at £ = 1.) Wahr et al. (1998) showed that a change is surface mass density A a(<j>, A) 
can be calculated from a change in the Stokes coefficients using 

= f ££(^Win<» 

e =2 m =0 ' £ ' 

x [AC* m cos mA + A S( m sin mA] ( 2 ) 

where p e is the average density of the Earth and k? is the elastic gravity load-Love number for degree 
£ (Farrell, 1972). 

If we assume that temporal changes the observed GRACE coefficients represent changes in the 
surface mass, then the gravitational load of this surface mass will deform the surface of the Earth. 
The elastic Green’s function G( 7 ) for radial displacement for a point an angular distance 7 from the 
point load (Farrell, 1972) is 
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where he is the elastic deformation load-Love number for degree £ and m e is the mass of the Earth. 
The radial deformation A r(<f>. A) for an arbitrary surface load A a(<j>, A) is 

A A) = a 2 JJdQ' A o{4>' , A') G( 7 ) (4) 

with the integration performed over the entire surface of the Earth, and where dfl' = dcfr'dX' cos <j>' 
and cos 7 = sin <p sin <j>' + cosc)cos</>'cos(A — A'). Using Eqs. (2) and (3) in Eq. (4), and using the 
Addition Theorem to express cos 7 , yields 
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3 Temporal variation of gravity 

Visual inspection of time series of observed Stokes coefficients leads us to conclude that there is 
a strong annual component to the temporal variation of gravity. With just slightly over a 2-year 
timespan, some apparently “secular” (or at least longer-term) variation is also visible in some of the 
time series. So that this longer-term variation does not systematically bias estimates of the annual 
variation, we include a rate term in our analysis. Thus, our models for the Stokes coefficients is 
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where Cf m is a constant value. Cf m is a rate term. C.^ m and are annua! amplitudes, and At = 
t - t a . (Superscript “i” stands for “in-phase” and “o” for “out-of-phase.”) The frequency / is one 
cy cle per year, and we used t 0 = 2003.0. 

Fig. 1 shows example time series for two Stokes coefficients, Cg,i and C 3 3.0- The error bars 
used are the formal errors from the Level-2 data files. The time series for Cg.i has a pronounced 
sinusoidal term as well as a visible secular variation, whereas the time series for C33.0 has neither. 

The fitting procedure enables us to evaluate the GRACE data noise. Following Wahr et al. (2004), 
we define the “realistic” errors SCe m and SS(m to be the weighted root-mean-square (vvrms) residual 
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The observed error degree amplitude is shown in the top frame of Fig. 2. The baseline GRACE 
values require a scaling of ~47 to achieve equality with the observed values. Wahr et al. (2004) 
found a value of 40 for the scaling, but they included degrees between 70 and 120, for w hich a 
smaller scaling is required (as can be seen from their figure). 

Using the scaling factor of 47, we can calculate realistic error degree amplitudes for the rate and 
annual-amplitude parameters of Eq. (6) and compare these to the degree amplitudes calculated using 
the parameter values. These are shown in the center (annual amplitude) and bottom (rate) frames 
of Fig. 2. For the seasonal terms, the estimated signal exceeds the error only for £ < 15 or so, in 
agreement with (Wahr et al.. 2004). The rate signal exceeds the error only modestly for all degrees 
except for £ = 2; the C9.0 term shows a large signal of currently unknown origin. 


4 Filtering the annual signal 

A common way of filtering the coefficients to enhance the contribution of a signal at low degrees 
and decreasing the contribution of noise at high degrees is to use an azimuthaliy sy mmetric Gaussian 
averaging function of Jekeli (1981). In the spherical harmonic domain, this function depends on 
degree only. However, we observe that the annual signal is not necessarily significant at all orders, 
even if the degree amplitude is significantly greater than the error. As an alternative, we explore a 
statistical approach. 

When the errors are not well known, we base the significance of the annual terms on a hy pothesis 
test for the significance of the improvement of the fit when the annual terms are included compared 
to when they are fixed to zero. If we have A 7 epochs, then the least-squares solution has N — 4 
degrees of freedom when the model Eq. (6) is used. If the annual terms are fixed to zero, the solution 
has N — 2 degrees of freedom. If < s the postfit y 2 statistic with N - 2 degrees of freedom and 
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X 2 with N — 4 degrees of freedom, we form the statistic 
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Under the hypothesis that the “true” model contains no annual terms, the numerator of the term in 
parentheses of Eq. (8) is the difference between x 2 statistics calculated for N — 2 and N —4 degrees 
of freedom, and is therefore has a x 2 distribution with 2 degrees of freedom. The denominator of 
the term in parentheses has a x 2 distribution with N — 4 degrees of freedom. Under the hypothesis, 
then, the $ has an F-distribution with (2, N — 4) degrees of freedom and an expectation of unity. 

A large difference between the x 2 statistic with and without the annual terms estimated would 
indicate that these terms “absorb” variation relative to the residual variation. This situation would 
lead to a large $ value compared to the expectation of unity under the hypothesis. If the annual 
terms absorb little or no variation, a small $ value would result. We can reject the hypothesis of 
no annual variations if the $ statistic is greater than some confidence limit for the F-statistic. We 
therefore choose only the set of (£, m) pairs for which we reject the null hypothesis. We will call 
this set fi, so that we use only £, m € 0. 

A problem presents itself in that we perform the fit and the F-test separately for the C and S 
Stokes coefficients. We have investigated two ways of handling this issue. In the first approach, we 
define the inclusion set Qc for the C coefficients and 0.s for the S coefficients. Another approach 
would be to define £, m £ Cl if the F-test for either the C or S coefficient indicated significance, 
yielding 
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Here the superscript a indicates either in-phase or out-of-phase annual amplitudes, calculated from 
our fit to Eq. (6). Because it does not make a significant difference for this data set, we will use this 
simpler approach for the following calculations. 


5 Observed annual deformations 

Our data consisted of 22 Level-2 GRACE data sets acquired between 2002.33 and 2004.54. Most of 
these data files contained Stokes coefficients with £ < 120, although several had only £ < 70. We 
therefore limited our analysis to £ < 70. Using the time series for each Stokes coefficient, we used 
least-squares to estimate the parameters of Eq. (6). We also performed the solution with the annual 
terms constrained to zero, so that the $ statistic of Eq. (8) could be calculated. We then formed 
the set fl of degree/order pairs that exceeded the 99.99% level of confidence for an F-distribution 
with (2,18) degrees of freedom. The $ statistic for 94 degree/order pairs exceeded this level of 
confidence. 
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Fig. 3 shows the observed global in-phase and out-of-phase annual radial deformation amplitudes. 
The largest deformation amplitudes are in areas where previous works (Tapley et al.. 2004; Wahr et 
al., 2004) indicate the largest amplitudes of inferred annual surface water amplitudes. The deforma- 
tion signal is largest in the Amazon Basin. sub-Saharan Africa, and southeast Asia. In the latter two 
regions, the amplitude of the deformation is ~6 mm. In the Amazon basin, the deformation signal 
is twice as large. 

The deformation field of Fig. 3 shows a number of smaller amplitude features, and the F-test 
filtering technique we applied does not yield for the largest features the round “bull's eye” shape 
that Gaussian smoothing yields (Tapley et al.. 2004: Wahr et al., 2004). In Fig 4 we show in larger 
scale the observed global deformations for South America calculated in three ways: the F-test filter 
describe in Sec. 4 (top) and Gaussian-smoothing of the complete set of coefficients with radii of 
1000 km (middle) and 300 km (bottom). (We will refer to the deformation amplitudes determined 
using the F-test procedure as “filtered” and those produced with Gaussian smoothing as “smoothed”) 
The Gaussian smoothing w ith a radius of 1000 km yields much "rounder" deformation patterns. The 
filtered deformations yields a deformation pattern that is highly elongated in the east-west direction. 
Because of this. Gaussian smoothing w ith a radius of 1000 km also yields smaller total deformation 
amplitudes, as the smoothing averages in values to the north and south that fall off to zero more 
rapidly. Gaussian smoothing with a radius of 300 km preserves the amplitude, but a great deal of 
“streaking” associated with GRACE errors is evident. The filtered reconstruction also shows 1-2 mm 
in-phase amplitudes throughout the southern portion of South America, whereas these amplitudes 
fall off in distance from the Amazon for the 1000-km-smoothed reconstruction. 

Davis et al. (2004) compared annual amplitudes of radial deformation from Global Positioning 
System (GPS) sites in South America to those determined from GRACE, and found a very high 
correlation. Here, w e repeat those calculations for the extended GRACE data set and also reconstruct 
the deformations using Gaussian smoothing. We then use the GPS-GRACE comparison to address 
the question of which of the three reconstructions (F-test filtering, Gaussian smoothing with 300 km 
radius, or Gaussian smoothing with 1000 km radius) yields more accurate amplitudes. 

To compare these two data types, we must first correct for the absence of the / = 1 term in the 
GRACE data. This correction procedure is described by Dav is et al. (2004). Briefly, we use the GPS- 
GRACE residuals to estimate three parameters representing the cartesian components of the annual 
amplitude of the geocenter variations, which express themselves as the f = 1 term. We then use these 
parameter estimates to “correct” the GRACE amplitudes for the £ = 1 term. The determination of 
the l = 1 corrections were done independently for the filtered and smoothed GRACE estimates. 

The GRACE- GPS comparisons are shown in Fig. 5. Each comparison displays a reasonably 
high correlation. The weighted root-mean-square GPS-GRACE difference is 1.5 mm for the filtered 
values, 1.8 mm for the 1000-km-smoothed, and 1.6 mm for the 300-km-smoothed. One large contri- 
bution to the better performance of the smaller smoothing radius compared to the larger smoothing 
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radius is that the larger amplitude seen in the GPS data is captured better. As discussed above, the 
larger smoothing radius averages out the larger amplitudes. This effect can be seen in Fig. 6, which 
shows the vvrms GRACE-GPS fit for a range of soothing radii. Starting at 1000 km, decreasing the 
smoothing radius improves the fit (decreases the vvrms) until, for radii less than ~300 km, the errors 
in the high-degree terms begin to dominate the comparison. However, visual inspection of Fig. 4 
indicates that a large amount of noise is apparent at a smoothing radius of 300 km. Thus, the fit of 
such a solution would depend critically on the geographical sampling represented by the GPS sites. 

The filter approach recovered the geographic pattern of annual deformation slightly better than the 
smoothing approach. It enabled the inclusion of higher degree terms that better define the shape of 
the deformation without introducing the “streaks.” However, a number of other larger-scale features 
appear (see Fig. 3) that we are unable to determine, without further investigation, whether they are 
real. 

6 Observed rates 

The time series of Cg,\ in Fig. 1 showed a definite linear trend. We might therefore ask whether any 
significant rates can be observed. Our model Eq. (6) contained a term representing a rate, and so it 
is a matter of reconstructing a rate field using these estimated parameters. The F-test method does 
not work reliably on this parameter, probably because there is only one degree of freedom (i.e., one 
parameter is being tested), and the rate estimate is extremely sensitive to outliers, especially near the 
beginning and end of such a short time series. We therefore choose to reconstruct the field using 
Gaussian smoothing, and we look at the geoid rate, since we expect to encounter secular variation of 
gravity not associated with elastic loading but also with glacial isostatic adjustment (G1A). However, 
we expect at this stage to see a number of other signals that are not truly secular but that are absorbed 
by a rate parameter linear over the two-year timespan of GRACE data. 

Although the observed geoid rates (Fig. 7) appear noisy, the observed values are well within 
reasonable bounds. In particular, the rate of change of the geoid over Hudson’s Bay in North America 
is ~1.1 mm yr -1 , consistent with the range of expectations (Wahr and Davis, 2002). The largest 
negative signal —1.6 mm yr -1 is over the Gulf of Alaska, where glacier melting is contributing 
~0.3 mm yr -1 of globally averaged sea-level change (Tamisiea et al., 2004). A longer time series 
will be required to determine the reliability of these and other features evident in Fig. 7. 

7 Conclusions 

The F-test “filter” performs better than Gaussian smoothing for the annual amplitudes in South 
America. It admits some high spatial frequencies into the solution and removes others whose coeffi- 
cients do not exhibit significant annual variation. As a result, for example, the east-west elongation 
of the deformation is observed that cannot be observed by smoothing. Preliminary studies indicate 
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that this approach works well on a global network for which the amplitudes are in general much less, 
although there are regional variations in correlation (Elosegui et al., 2004). 

Although the rates estimates, which result from the same solutions as the annual amplitudes, are 
noisy, at least one known signal can be observed. GIA in Hudson's Bay Of course, the goal of 
GRACE is to measure unknown signals. It is difficult to know at this early stage, however, which of 
the signals are truly secular. 
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Figure Captions 


Fig. 1. Example time series for GRACE Stokes coefficients (shown as ±1 <t error bars) and best-fit to Eq. (6) 
(gray line). 


Fig. 2. Top: Degree amplitudes for the GRACE errors. The blue line are the observed degree amplitudes, 
the dotted red line are the GRACE baseline values, and solid red line is the GRACE baseline scaled by 47. 
Middle: Degree amplitudes for the errors in the annual amplitudes (red) compared to the degree amplitudes for 
the estimated parameters (blue). Bottom: Degree amplitudes for the errors in the rate (red) compared to the 
degree amplitudes for the estimated parameters (blue). 


Fig. 3. Observed in-phase (top) and out-of-phase (bottom) annual radial deformation amplitudes, reconstructed 
using the observed annual amplitudes and Eq. (9). 


Fig. 4. Larger-scale figure focusing on annual amplitudes for South America. The top row is calculated using 
the F-test filtering of Eq. (9) The middle row is calculated using Gaussian smoothing with a radius at half- 
maximum of 1000 km, and the bottom row used a radius of 300 km . The left column is in-phase, and the right 
column is out-of-phase. Also shown are the locations of the GPS sites used in the comparison. The color scale 
is the same as that of Fig. 3. 


Fig. 5. Comparison of GPS and GRACE-determined annual amplitudes using the F-test filter (top) and Gaus- 
sian smoothing with radius 1000 km (middle) and 300 km (bottom). The GRACE amplitudes have been cor- 
rected for the absence of the ( = 1 term in the GRACE data as described in Davis et al. (2004), which also 
describes the analysis of the GPS data. 


Fig. 6. Weighted rms (wrms) fit for GPS-GRACE comparison for a range of smoothing radii. The filtered 
amplitudes yielded a wrms fit of 1 .48 mm 


Fig. 7. Observ ed geoid rates using Gaussian smoothing with 1000 km radius. 
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Fig. 2. Top: Degree amplitudes for the GRACE errors. The blue line are the observed degree amplitudes, 
the dotted red line are the GRACE baseline values, and solid red line is the GRACE baseline scaled by 47. 
Middle: Degree amplitudes for the errors in the annual amplitudes (red) compared to the degree amplitudes for 
the estimated parameters (blue). Bottom: Degree amplitudes for the errors in the rate (red) compared to the 
degree amplitudes for the estimated parameters (blue). 
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Fig. 3. Observed in-phase (top) and out-of-phase (bottom) annual radial deformation amplitudes, reconstructed 
using the observed annual amplitudes and Eq. (9). 
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Fig. 4. Larger-scale figure focusing on annual amplitudes for South America. The top row is calculated using 
the F-test filtering of Eq. (9). The middle row is calculated using Gaussian smoothing with a radius at half- 
maximum of 1000 km. and the bottom row used a radius of 300 km . The left column is in-phase, and the right 
column is out-of-phase. Also show n are the locations of the GPS sites used in the comparison. The color scale 
is the same as that of Fig. 3. 
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Fig. 5. Comparison of GPS- and GRACE-determined annual amplitudes using the F-test filter (top) and Gaus- 
sian smoothing with radius 1000 km (middle) and 300 km (bottom). The GRACE amplitudes have been cor- 
rected for the absence of the i — 1 term in the GRACE data as described in Davis et al. (2004), which also 
describes the analysis of the GPS data. 
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Fig. 6. Weighted mis (wrms) fit for GPS-GRACE comparison for a range of smoothing radii. The filtered 
amplitudes yielded a wrms fit of 1 .48 mm 
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Fig. 7. Observed geoid rates using Gaussian smoothing w ith 1000 km radius. 
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